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Delayed outflows from black hole accretion tori following 
neutron star binary coalescence 
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ABSTRACT 

Expulsion of neutron-rich matter following the merger of neutron star (NS) binaries is 
crucial to the radioactively-powered electromagnetic counterparts of these events and 
to their relevance as sources of r-process nucleosynthesis. Numerical simulations of 
NS-NS coalescence find, however, a wide range in the quantity of prompt dynamically- 
ejected mass. Here we explore the long-term (viscous) evolution of remnant black hole 
accretion disks formed in such mergers by means of two-dimensional, time-dependent 
hydrodynamical simulations. The evolution of the electron fraction due to charged- 
current weak interactions is included, and neutrino self-irradiation is modeled as a 
lightbulb that accounts for the disk geometry and moderate optical depth effects. 
Over several viscous times (~ 1 s), a fraction ~ 10% of the initial disk mass is ejected 
as a moderately neutron-rich wind (Y e ~ 0.2) powered by viscous heating and nuclear 
recombination, with neutrino self-irradiation playing a sub-dominant role. Although 
the properties of the outflow vary in time and direction, their mean values in the heavy- 
element production region are relatively robust to variations in the initial conditions 
of the disk and the magnitude of its viscosity. The outflow is sufficiently neutron-rich 
that most of the ejecta forms heavy r-process elements with mass number A > 130, 
thus representing a new astrophysical source of r-process nucleosynthesis, distinct 
from that produced in the dynamical ejecta. Due to its moderately high entropy, disk 
outflows contain a small residual fraction ~ 1% of helium, which could produce a 
unique spectroscopic signature. 

Key words: accretion, accretion disks — dense matter — gravitational waves - 
hydrodynamics — neutrinos — nuclear reactions, nucleosynthesis, abundances 
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Coalescing neutron star (NS) and stellar mass black hole 
(BH) binaries are among the most promising sources for de- 
tection with networks of ground-based gravitational wave 
( GW) observatories, inc l uding Advanced LIGO and Virgo 
Abramovici et al.l 11992; lAcernese et al.l |200S| ; ICaron et al.l 



Il999l ; lAbadie et all l20ld . l2012l l. The GW signal pro- 
vides infor mation on the p r operties of the m erging bi- 
nary (e.g., lAiith et all 120071 ; Ivan der Sluvs et al.l [2008 1; a 
potentia l probe of the equation of state of dense ^mat- 
ter (e.g.lFaber et al.ll2002l ; Ishibatal 120051 ; iRead et all 120091 ; 
iBauswein et al.l 201211; and eve n a test of General Relativity 



itself fe.g. lCornish et ai]l201ll l. 

To maximize the scientific benefit of a GW detec- 
tion, it is important to identify a coincident electromag- 



netic counterpart (e.g . Bloom et al.ll2009l . Metzger fc Bergerl 



|2012| . INissanke et allfeoil iBartos et al.ll2012l 1.The sky lo- 



calization provided by initial Advanced L IGO /Virgo will 
be tens to hundreds of square degrees (e.g. iFairhurstl 120091 ; 



I Wen fc Chen1l2010l ; INissanke et al"]|201ll 1. but a coincident 
X-ray, optical or radio signal could provide a much more 
accurate position, such that the host galaxy and redshift 
of the merger could be identified. Electromagnetic emission 
also indirectly probes the (magneto-)hydrodynamics of the 
merger and its aftermath. 

The BH created from a NS-NS merger, following a 
metastable hypermassive NS phase, is typi cally surrounded 
by a t hick torus of mass ~ 10~ 3 — O.IMm IjRuffert fc J anka 
19991; lUrvuetall I2OO0I; iRosswog fc Liebendorferl 120031 



Oechslin fc JankalhoOrj; 
2010; iHotokezaka et al.l 



Ghawla et al.ll20ld ; iRezzolla et all 
forms follow- 



2011 



. A torus 
binary 



also 



ing a BH-NS merger if the binary mass ratio is suffi- 
ciently low that the NS is tidally disrupt ed outside the BH 
instead of being swallowed whole (e .g. Faber et all 2006; 



Duezetal. 201 



_ . _ _ . _.. .Foucart et alj|201l] ; Istephens et al 

Foucart et al.ll2012l ; lEast et al.ll2012l l 



2011; 



The short accretion 
and large accretion rates M > Mq 



timescale t v isc < 1 s, 

S~ , of such remnant tori motivate the current paradigm 
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that NS-NS/NS-BH mergers are the central engi nes power- 
ing short duration gamm a-ray bursts (SGRBs) (jPaczvnskil 
ll986l : lEichler et al.lll989D . 

Despite being well-studied, SGRBs do not necessar- 
ily represent an ideal counterpart for most GW-detected 
mergers because SGRBs are difficult to localize on the 
sky and their observed rate is relatively low (< 1 yr _1 ) 
within the sensitivity r ange of Advanced LIGO/Virgo 
JMetzger fc Bergerll20lj ; IChen fc Ho"izll2012l ). These draw- 
backs have motivated the study of more isotropic merger 
counterparts, such as supernova-like optical transients pow- 
ered by the decay of radioa ctive elements synthes i zed in the 
merg e r ejecta ('kilonovae' , iLi fc Paczvhskil 1 19981 ; iKulkarnil 
120051 ; IMetzger et all l20ld ). Merger simulations show that 
a modest quantity of highly neutron rich material (electron 
fraction Y e < 0.1) is ejected by tidal forces during the merger 
(ILattimer fc Schramm! I1974: Rosswog et al J 11999: iRosswoe 



Chawla et al.l |2010l: IStephens et al.ll201ll: lEast et al 



2012 : iRosswoa 120121 : iRosswog et alj 120131 : iBauswein et al 
2013J), and recent nuclear reaction network calculations 



largely agree on the amount of radioac tive heating provided 
by such decaying r-process elements ( Metzger et al.l l20ld : 
iRoberts et al]l201ll : lGorielv et al.l[201ll ). 

One large uncertainty, however, is the opacity of the 
ejecta, because current kilonova models predict a composi- 
tion dominated by heavy r-process nuclei (atomic mass num- 
ber A > 130). Recent atomic structure calculations show 
that the opacity of material containing even a small quan- 
tity of Lathanide elements (A ~ 139 — 175) could be or- 
ders of magnitude hig her than that of Fe (jKasen et alj|2013l ; 
iBarnes fc Kasen 2013). This implies that the predicted tran- 
sient peaks on longer timescales - up to a ~ week - and at a 
lower lumino sity and frequency - in the near infrared rather 
than optical (jBarnes fc Kasenll2013l ) - than i n original mod- 
els that assumed opaci ty similar to Fe (e.g. iMetzger et al.1 
|2010| . IPiran et al.ll2013h . Given the sensitive dependence of 
the signal on the Lathanide fraction, it is thus crucial to 
identify any diversity in the nucleosynthetic composition of 
the ejecta. 

Another uncertainty is the mass of the ejecta M e j. 
Although the quantity of matter expelled in dynamical 
'tidal tails' from the merger can be l arge (> W~ 2 Mc.) ) 
in some cases (e.g. eccentric mergers; lEast et al.l I2012T ). 
it is ofte n found to be significa ntly smaller (~ 10 -4 — 
10~ 3 Mm; iHotokezaka et aljBoill ). Some fully general rel- 
ativistic merger simulatio ns find that essentially no mass at 
all is ejected dynamically (jFaber fc Rasidl2012l ). 

Dynamical expulsion is not the only source of ejecta, 
however. Mass loss also occurs in outflows from the accre- 
tion disk over longer, viscous timescales. Initially, neutri- 
nos radiated from the disk - or from the central hyper- 
massive NS prior to its collapse - can heat and drive an 
outflow from the surface of the disk (a 'neutrino-driven' 



wind; e.g. McLaughlin & Surman 2005; Metzger et al. 
Surman et al. 2008; Dcs sart et al.M2009l : iKizivat et al. 



2008; 



2010; 



Wanaio & Janka 2012J). A potentially even larger quantity 
of mass is lost from the disk at later times due to power- 
ful outflows driven by viscous heating and recombination 
of free nucle i into q-particles (ILee fc R amirez-Ruiz 2007; 
Beloborodovl 120081; IMetzger et all |200S| , l2009al : iLee etlil 



2009J ) 



a height-integrated, time-dependent model of the viscous 
spreading of accretion disks formed from NS-NS/NS-BH 
mergers, including the evolution of the midplane electron 
fraction Y e due to weak interactions. MPQ09 find that as the 
disk spreads due to angular momentum transport and the 
temperature decreases, e~/e + captures become slow com- 
pared to the viscous evolution timescale. Since these cap- 
tures both cool the disk and affect its composition, the disk 
becomes geometrically thick and Y e freezes out. Soon after 
freeze-out, free neutrons and protons recombine into a par- 
ticles, starting near the outer edge of the disk and moving 
inwards with time. MPQ09 estimated that ~ 20 - 50% of 
the initial disk mass is unbound by this energy deposition, 
with a range of electron fractions Y e ~ 0.1 — 0.4 result- 
ing from varying local conditions during freez e-out (see also 
iLee fc Ramirez-Ruiz|[2007l:lBeloborodovil2008l'). 



Lee et al.1 1 20091 ) performed axisymmetric (2D) global 



IMetzger et~all |2009al . hereafter MPQ09) constructed 



simulations of long term disk evolution, which confirmed 
this process of He 'evaporation' explicitly. However, since 
they assumed /3-equilibrium instead of explicitly following 
the weak interactions, they could not accurately determine 
the energy released by a-recombination or the electron frac- 
tion of the final ejecta. The electron fraction of disk out- 
flows is critically important because it controls whether the 
heavy elements synthesized in the outflow are dominated by 
2nd/3rd peak r-process elements (if Y e < 0.3 — 0.4); lighter 
neutron-rich elements (if 0.3 — 0.4 < Y e < 0.495); or sig- 
nifica nt quantities of 56 Ni (Y e > 0.495) (e.g. iHoffman et~al1 
ll997j). Given the potentially large difference in the optical 
opacity of Fe-group and r-process elements, and the differ- 
ent geometry of tidal tails and disk winds, the resulting kilo- 
nova light curve, including contributions from both sources 
of ejecta, could in pri nciple be much more complex than 
previously anticipated (jBarnes fc Kasenll2013h . 

Beyond their implication for electromagnetic emis- 
sion, disk outflows from binary NS mergers also rep- 
resent an essentially unexplored astrophysical site for 
heavy element synthesis. This issue is particularly im- 
portant because the astrop hysical origin of r - proces s 
nuclei remains a mystery IjQian fc Wasserburd 120071 ). 
and binary NS mergers c ould be an importan t — and 
possibly dominant— sou r ce (ILattimer fc Schramm! |l974j; 
iFreiburghaus et al.1 Il999l : iKorobkin et al.ll2012l ). Neverthe- 
less, with the exception of neutrino-driven disk outflows 
(e.g. ISurman et ail J2008l : IWanaio fc Jankal |2012| ). to date 
most studies of the r-process in such environments focus 
on the dynamical ejecta. 

In this series of papers, we explore the long term evolu- 
tion of the accretion disks formed by NS-NS/NS-BH mergers 
by means of axisymmetric (2D) time-dependent hydrody- 
namic simulations. In paper I (this work) we explore the 
dynamics of the disk evolution and the properties of the re- 
sulting disk outflows. In separate works, our results for the 
disk outflows will be used to make predictions for the result- 
ing transient electromagnetic signature and nucleosynthetic 
yields. 

Our models employ a finite-volume hydrodynamic 
method, a realistic equation of state, and a complete im- 
plementation of charged-current weak interactions in the 
optically-thin regime. Approximations are made in order 
to capture the key physical processes at a low computa- 
tional cost: (1) the effects of general relativity on the disk 
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dynamics are approximated via a pseudo-Newtonian poten- 
tial, (2) neutrino self-irradiation is modeled as a 'lightbulb' 
that accounts for the basic disk geometry and includes cor- 
rections for moderate optical depth effects, and (3) angu- 
lar momentum transport is approximated by an anomalous 
shear-stress tensor with an a-viscosity. These approxima- 
tions enable us to evolve accretion disks for a large number 
of orbits (~ 1000) and thus perform a parameter space study 
on the outflow properties which can guide future, more de- 
tailed calculations. 

The paper is organized as follows. In S|2] we briefly sum- 
marize the properties of the initial accretion disk. The nu- 
merical implementation is described in S|3] Results are pre- 
sented in ij4l followed by a discussion in S|5] Our conclusions 
are summarized in SjB] The appendices provide additional 
information on the neutrino treatment as well as general in- 
formation on accretion regimes and nucleosynthesis in NS 
merger disks. 



2 DISK PROPERTIES 

The mass of the torus formed during a NS-NS merger can 
vary considerably, M t ~ 10 -3 — 0.1M©, depending on the 
binary mass ratio, the assumed EOS, the eccentricity of 
the initial orbit, and whether the system undergoes a long- 
lived hypermassive NS phase prior to BH formation (see 
iFaber fc Rasioll2012i for a recent review). The disk mass in 
a NS-BH merger also depends on the spin of the black hole, 
since spin controls the location of the innermost circular or- 
bit. If the NS is tidally disrupted outside the BH horizon, 
then the torus mass is relatively large M t ~ O.IMq; other- 
wise, the NS is swa llowed whole and little or no torus forms 
fe.g. . lFoucartll2012r ) . Merger simulations find that the initial 
torus is distributed across a wide range of radii, with most of 
the mass and angular momentum concentrated on a radial 
scale R ~ 20 - 100 km. 

The (Newtonian) orbital time at this radial scale is 



t„ 



3.5 x 10" 



R 







50 km 



3/ 2 



/Mb 



\3M e 



-1/2 



(1) 



where Mbh is the BH mass. The characteristic timescale for 
matter to accrete may be estimated by the viscous time 

1/2 
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(2) 



where Ho is the vertical scaleheight and a parametrizes the 
disk viscosity f £|3.3p . resulting in a characteristic accretion 
rate 

M t 

£visc 



Ah 



10~^M Q s 



V0.03/ 



M t 



0.03/ V 10 ~ 2M © 



R 







50 km 



-3/2 
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(3) 



V3M , 
Due to the high densities and rapid evolution of 



the disk, photons are not an important source of cool- 
ing on timescales of relevance. However, the tempera- 
ture is sufficiently high that neutrinos are relevant, espe- 
cially in the earliest p hases of the disk evolution (e.g.. 



Popham et~aill 19991; iNaravan et all l200ll ; |Pi Matteo et al I 



20021 : IChen fc Beloborodovl 120071 : hereafter CB07). A 



straightforward calculation shows that a thin neutrino- 
cooled disk (H ~ 0.1 — 0.3i?o) obtains when the accretion 
rate at radius Ro ex ceeds a critical value M„ given by (see 
iMetzger et alll2008l , their eq. [11]) 



M, w 1.3xlO _a Af Q B _1 (J^- 
\50km 



5/6 



a N 5 / 3 
O03 



Mb 



3M G 



1/2 



.(4) 



where the dominant source of cooling is assumed to be 
charged-current weak interactions, and general relativistic 
corrections have been neglected. Unlike photon-cooled disks, 
such as those in X-ray binaries or AGN, neutrino-cooled 
disks are the most radiatively efficient at small radii. 

By equating Ro in equation (|4|) with the radius i?i SC o of 
the innermost stable circular orbit of the BH, one obtains a 
critical accretion rate Mi gn (the 'ignition' rate) above which 
the disk is thin near its innermost radii: 



Mu 



iXl0 ^ 2xl0 "Hcl03 



5/3 /AfeH^ 4/3 

V3M e 



M s-\(5) 



where the prefactor was calculated by CB07 for a 3M© black 
hole of spin a = 0[a = 0.95] using a steady-state accretion 
model in a Kerr metric. 

When the disk is neutrino-cooled (M > M v ) its in- 
ternal energy eint is dominated by non-degenerate nucle- 
ons (CB07). The characteristic temperature of the torus 
midplane, which is attained on a cooling timescale t C ool ~ 
t v isc(H/Ro) , is thus approximately given by 



,,.-... + ):.(£)(£Vr!5-Wv. 



V Ro ) 



(6) 



where Y e — n p /(n n +n p ) is the electron fraction and n p (n n ) 
are the number densities of protons(neutrons) , respectively. 
The initial temperature of the disk is sufficiently high that all 
nuclei are dissociated into free nucleons. Alpha particles and 
heavier nuclei form only once the disk temperature decreases 
to < 1 MeV. 



3 NUMERICAL SETUP 

3.1 Equations Solved and Numerical Method 

We solve the equations of mass, momentum, energy, and 
lepton number conservation on an axisymmetric spherical 
polar grid (r,9), with source terms due to gravity, shear 
viscosity, and optically-thin weak interactions 



|+V-(pv p ) 







*£ = t-ivp-v* 

at p 

At 

p —± = rsin6»(V-T) 

p '" + pV • v„ = — T : T + pQnot 
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dt 



(7) 
(8) 
(9) 

(10) 

(11) 
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Here p, p, v p = v r f + veO, and eint denote the fluid density, 
pressure, poloidal velocity, and internal energy, respectively. 
The Lagrangian differential operator is d/dt = d/dt + v p - V. 
The specific angular momentum along the symmetry axis is 
given by l z — rsinOv^, with V4, the azimuthal velocity, and 
f c is the centrifugal force in the poloidal direction. 

The pseudo-n e wton ian gravitational potential of 
iPaczvhskv fc Wiital l|l980h is employed, $ = -GM B n/(r - 
Rs), with Rs = 2GMbh/c 2 the Schwarzschild radius, as ap- 
plicable to a BH of moderate spin[j The black hole mass is 
kept constant in the simulation. The neutrino source terms 
are Q not and r nct (see £|3.2p , and the viscous stress tensor is 
denoted by T (see H3.3P ^_^ 

The Helmholtz EOS (JTimmes fc Swesty|200ol) is used to 
relate internal energy and pressure. The ion component is an 
ideal gas of neutrons, protons, and alpha particles, with rel- 
ative abundances that satisfy nuclear statistical equilibrium 
(NSE). Heavier nuclei are not included, since the energy re- 
lease in their formation is much smaller in comparison, and 
hence should have little impact on the outflow dynamic^. 
An additional (inert) hydrogen gas is used to populate the 
ambient medium that initially surrounds the torus. The zero 
point of energy is taken to be the pure nucleon state, with 
the internal energy becoming 



eint — eint.O An 



(12) 



where ei n t,o is the specific internal energy provided by the 
Helmholtz EOS, X a is the mass fraction of alpha particles, 
m a is the mass of an alpha particle, and Q a — 28.3 MeV 
is the nuclear b inding energy of an alpha particle (e.g., 
I Audi et al.l 120031 ). In addition, the thermodynamic deriva- 
tives of the pressure, internal energy, and entropy of the 
ions acquire terms that originate in the dependence of the 
ion abundances on temperature and density. If the temper- 
ature falls below 5 x 10 9 K, or if the mass fraction of inert 
hydrogen exceeds 1%, the abundances are frozen. Any hy- 
drogen mixed with torus material is added to the proton 
mass fraction if the temper ature is highe r than 5 x 10 9 K. 

We use FLASH3.2 (iDubev et all I2009T ) to evolve 
the system of equations 0-JHJ with the dimensionally- 
split version of the Piece wise Parabolic Method (PPM, 
IColella fc Woodward! 1 1984 ) . The public version of the code 
has been modified to allow for a non-uniformly spaced grid 
in sph erical polar coordinates, as described in iFernandej 
d2012h . 

The computational grid is logarithmically spaced in ra- 
dius, covering the range 2Rs to 2 x 10 3 i?s in most models. 
In the polar direction, the cell spacing is uniform in cos8, 
covering the interval [0, n]. We adopt this grid structure be- 
cause we intend to cover a large dynamic range in radius, 
and to concentrate the resolution near the midplane, given 
the quasi-spherical character of the outflows. Our standard 
resolution uses 64 cells per decade in radius and 56 cells in 
the polar direction, yielding approximately square cells at 
the midplane (Ar/r ~ A6 ~ 2°). The coarsest angular cell 



1 BHs formed from NS-NS mergers ha ve typical spin parameter 
a - 0.7 - 0.8 Ce.g.. lRezzolla et aljlioich 

2 The formation of heavy nuclei would indeed occur at radii 
within our simulation grid, however, see i|5.1l 



next to the polar axis is 10.8°. One model is evolved at dou- 
ble the resolution in both radius and in angle, to quantify 
convergence. 

The boundary conditions at the polar axis are reflecting 
in all variables. At the inner and outer radial boundaries, we 
impose a zero-gradient boundary condition for all variables. 
In addition, the radial velocity in the ghost cells of the ra- 
dial boundaries is set to zero if its value in the first active 
cell next to the boundary would imply mass entering the 
domain. 

All energy and viscous source terms are set to zero when 
the density falls below a fiducial value of 10 g cm - . This 
value is sufficiently low as to not impact the outflows from 
the torus, yet sufficiently high to avoid requiring too small of 
a time step due to angular momentum transport in (unim- 
portant) low-density regions. 

To prevent numerical problems in the funnel that de- 
velops around the rotation axis, we impose a floor of density 
with a value ~ 20% lower than the ambient density, 1 g cm -3 
(j HH) , ano - a floor of temperature at 10 K. 



3.2 Neutrino Treatment 

The inner regions of the disk become opaque to neutrinos 
above an accretion rate (Appendix |A")| 



M 



O.lSMes- 1 (_2L^ 6 (J£ 

V0.03/ \3M, 



7/10 



R 



9/5 
'-"' '' .(13) 



6R 



The condition M 



> Mo (eq. [3]), sets the minimum 
disk mass for optical thickness. For slowly spinning BHs, 
this mass is ~ O.IMq. 

Most of the disks we explore (Mt < O.IMq) are thus 
only marginally optically thick to neutrinos at early times. 
Since the disk density decreases as accretion proceeds, full 
neutrino transparency is always achieved. The models we 
explore undergo this transition before weak freezout occurs, 
as quantified by the magnitude of the optical depth, and 
hence we adopt an optically thin treatment of energy ex- 
change with matter and regulation of the electron fraction 
in all cases. Finite optical depth corrections are included to 
prevent excessive heating at early times in the few massive 
tori we evolve. 

Charge d-curren t weak interactions are inclu ded as for- 
mulat ed by iBruennl (|l985T ) and implemented by IFernandej 
1)20121) . These reactions control the evolution of the electron 
fraction, and provide the dominant energy exchange chan- 
nel between neutrinos and matter. The relevant source terms 
are the net rate of change of the electron fraction 



1 net — J- e — J- e 



(14) 



where F e _ and F e + are the rates of electron and positron 
creation per baryon, respectively, and the net specific heat- 
ing rate of the gas 



— iLi/p 1 iLf e ^si/p \~s~Vp ^P 



(15) 



where H Ve and Ha c is the heating by absorption of neutri- 
nos and antineutrinos on nucleons, respectively, and C„ e and 
Cj> e are the cooling rates from electron and positron capture, 
respectively. Explicit forms for these terms are provided in 
Appendix [B] Additional neutrino energy losses from pair 
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creation (C pa irs) are included via the lltoh et al.l (|1996f ) pa- 
rameterization^ Source terms are included via explicit op- 
erator split in the equations for energy and lepton number 
conservation. 

Self-irradiation of the disk is implemented via a 
lightbulb-type prescription that accounts for the disk geom- 
etry (Appendix |BJ. All of the neutrino emission is assumed 
to come from a ring of material on the midplane, with a 
radius -Rem chosen to be an average value weighted by the 
neutrino emissivity, 

fsm9d9drr 3 p{C„ c +CpJ 



R 



Jsm9d9drr 2 p(C v<! +Cp e )' 



(16) 



This approximation is justified by the highly peaked form 
of this emissivity (Appendix [B} , and the low optical depth 
of most of the disk masses considered. The neutrino lumi- 
nosities of electron-type neutrinos L„ c and antineutrinos L^ c 
used in this self-irradiation treatment are calculated at the 
end of the previous time step by integrating the correspond- 
ing emissivities over the whole torus. Spectra are assumed 
to follow a Fermi-Dirac distribution with zero chemical po- 
tential, with a temperature equal to the gas temperature 
in an angle-averaged radial shell (10% spread) around the 
emission radius Rem- Self-irradiation is included to explore 
its influence on the ejecta composition and its relative en- 
ergetic contribution to the total mass loss. Results are not 
expected to be very sensitive to the particular choices of lu- 
minosities and spe ctra, because the proce ss is energetically 
subdominant fe.g-. lRuffert fe Jankalll999l . MPQ09). 

To explore more massive tori, finite optical de pth effects 
are ap proximated by a modified version of the I Lee et al.l 
i|2004r ) treatment. At each time step, four optical depth func- 
tions for each neutrino species are computed from the neu- 
trino emission peak, 

r e 

t? p (0) = / K t R cm d9', 9e [0,71-/2], r = R am (17) 

Jo 



tT(0) 



out / \ 

n (r) 



Ki flemdfl', 9 e [tt/2, 71"], T = Rem (18) 



Kidr', r e [Rem , r out ] , 9 = — (19) 



Ti n {r) = 



Ki dr', 



re [n„, Rem], 9=~, (20) 



corresponding to the four coordinate directions. Given that 
most of the torus mass lies within ~ 45° of the equator, 
we use angular integration to estimate the vertical optical 
depth. For fast computation, the absorption coefficie nts for 
this o ptical depth functionqj are approximated as in Ijankal 
(|200lT ) 



I" 1 l'. 



2.2 x IQ' 7 p w T^ A y/X p {2X n +X P ) cm" 1 , (21) 
2.2 x 10~ 7 p w T^ A ^/X v {2X v +X n ) cm" 1 . (22) 



The neutrino emissivities Cu c and &> e are then suppressed 
everywhere by a factor exp (— Tj, pea k), where 

Ti.pcak =mhl[T^ P (Tr/2),T? n (n/2),T° nt (Rem),Tr(Rem)]- (23) 



3 http://cococubed.asu.edu/codc_pagcs/codcs.shtml 

4 Not to be confused with the absorption coefficients used in the 
actual neutrino absorption rates (Appendix FBI) , 
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This prescription accounts for the fact that most neutri- 
nos will escape along the path with the least optical depth, 
which in most c ases is the vertical direction. In contrast to 
iLee et al.l (|2004r i. we ignore the neutrino pressure contribu- 
tion, which is a small correction (our tori are not radiation- 
dominated while optically thick). 

The neutrino luminosity used for self-irradiation, first 
computed from the attenuated emissivities, is further sup- 
pressed by a local factor exp(— Ti,m), where 



Ti,irr(r, 9) 



= [rf (-Rem) - T, 

, T?(r) 

T?(Ren 



(r)] 
f(n/2)-rf(9)}. 



(24) 



The superscripts x and y are the radial and angular optical 
depths for the appropriate quadrant, respectively (eqns [17] - 
[20]). This prescription attempts to capture the attenua- 
tion of the neutrino flux in regions far from the emission 
peak, and smoothly approaches the semi-transparent light- 
bulb limit for decreasing torus density. It is a coarse approx- 
imation, however, since it does not account for attenuation 
along lateral directions, and assumes that the vertical den- 
sity and abundance dependence differs only by a constant 
scaling factor at different radii. 

In optically thick disks, the rate of change of the elec- 
tron fraction (eq. [14]) is affected indirectly via the suppres- 
sion of the self-irradiation luminosities. The portion due to 
the local neutrino emissivities is not modified, however, since 
its contribution will still vanish in beta equilibrium. 



3.3 Angular Momentum Transport 

Angular momentum transport is mediated by a viscous 
stress tensor T fo r which only the azimuthal components 
are non-zero fe.g-. IStone et al.lll999n . 



T r <j> — pV 



Tetf, — 



pv 



r d 


sin 9 dr 


sin 9 d 



(25) 
(26) 



Including all components of the stress would sup- 
press convection in the poloid al direction (e.g., 
Ilgumenshchev fc Abramowica Il999f ). Our simulations 
mimic turbulent angular momentum transport via 
thermally-driven convection. Given the low mass of 
the torus relative to that of the central BH, and its 
relatively thick scaleheight, we consider our neglect of self 
gravity as a good approximation. 

To connect wit h prev ious calculations, we employ a 
IShakura fc Sunvaevl |l973j) parameterization of the kine- 
matic viscosity coefficient 



Q f 



GM BH /r 



VGMbhRo 1 



R* 



Ro 



1/2 

(27) 



where fit 



s t he K eplerian frequency of the 
IPaczviiskv fc Wiital |l980j) potential and cf = p/p is 
the isothermal sound speed. 

The implementation of angular momentum e v olutio n 
in FLASH is described in iFernandez fc Metzgeil l|2013h . 
Viscous heating is handled as an explicit source term in 
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Table 1. Models Evolved and Outflow Properties. Columns from left to right are: initial torus mass, black hole mass, initial torus radius 
(eq. 1281 ), initial electron fraction, initial entropy, viscosity coefficient (eq. 1271 ), spatial resolution (std=standard, hi=high; £|3. It . inclusion 
of recombination energy (eq. 1121 ), inclusion of self-irradiation f i|3,2l l, mass-flux- weighted values of the electron fraction, entropy, and 
expansion time t = r/v r at a radius where f ~ 5 X 10 9 K and within ±60° of the midplane (eq. 1311 ), ratio of the total mass with 
positive energy ejected beyond r = 10 9 cm to the total mass accreted through the ISCO radius, each in units of the initial torus mass, 
and mass-flux-weighted density (in units of 100 g cm -3 ) and velocity (in units of 10, 000 km s — 1 ) at r = 10 9 cm and within ±60° of the 
midplane. 
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the energy equation, while the diffusive momentum flux is 
added to the advective fluxes obtained during the Riemann 
solver step. Tests of t his im plementation are provided in 
iFernandez fc Metzger! J20131 ). 



3.4 Initial Conditions 

The initial condition is an equilibrium torus with constant 
entropy, angular momentum, and electro n fraction (e.g., 
iPapaloizou fe Pringld Il984l ; iHawlevI I2OO0T ). The initial so- 
lution is obtained by first finding the local density given the 
entropy, electron fraction, and position in the disk, 



w(p,r,8) 



GMbh 
Ro 



Ro 



Rs 



ii'n 



2(rsinr?) 2 



Ro 



R, 



2d 



(28) 



where w = a n t + p/p is the specific enthalpy of the fluid, 
and d is the torus distortion parameter. The resulting den- 
sity distribution is then integrated to obtain the disk mass. 
The procedure is repeated, changing the entropy, until the 
desired torus mass is obtained. The distortion parameter d 
controls the thermal content of the solution; for fixed mass 
and electron fraction, it is uniquely determined by demand- 
ing a given entropy or a value of H/R = Ci/(rf2x) at the 
point of maximum density. 

Evolving this equilibrium torus against a low-density 
background results in the emergence a shock from the edges 
of the disk due to strong gradients. To avoid confusing this 
initial transient with the viscously-driven outflow we are in- 
terested in, we evolve the initial torus for 100 orbital times 
at r = Ro without angular momentum transport or energy 
source terms. At this time a power-law density profile has 
been established at large radii, and the total mass in torus 



material has decreased by a small amount (< 1%), due to 
to the spread to larger radii. The mass accretion rate at the 
ISCO due to numerical viscosity is ~ 10 _5 Mq s _1 . We then 
reset all positive radial velocities to zero and take this as the 
initial condition for the simulations. 

A realistic merger will produce a torus with a more 
general radial distribution of mass, angular momentum, and 
thermal energy, than we have assumed. However, we do not 
expect our results to be overly sensitive to the initial thermal 
distribution of the torus, since the temperature profile will 
equilibrate on a thermal timescale t coo i ~ t v i BC (H/R) 2 ~ 

U. l£visc - 

The torus is surrounded by an ambient medium of con- 
stant density equal to p a mb = 1 K cm ~ 3 i an d pressure distri- 
bution p am b = — $Pamb l|Stone et al.lll999l ). The magnitude 
of the ambient density is chosen as the lowest value that 
does not cause numerical problems in the inner regions of 
the polar funnel. 



3.5 Models Evolved 

We evolve a number of models that explore the effect of (1) 
varying fundamental parameters in the system, and (2) ar- 
tificially removing critical processes for outflow generation. 
All models are shown in Table [l] 

The S-series of models include all the physics, and vary 
parameters around a fiducial model with Mbh = 3Mq , a = 
0.03, M t - 0.03M Q , Y e = 0.1, and s = 8 k s per baryon, with 
standard resolution ( i]3.1|) . We call this model S-def, and it 
is intend ed to roughly ma t ch th e parameters of model B10 
PaWi of lRuffert fc Jankal (|l999l . cf. their Figure 14). This 
model is also marginally optically thick to neutrinos at early 
times, becoming transparent (as measured by the minimum 
optical depth r p0 ak [eq.[23]) after 5 orbits at r — Ro. Models 
in this series explore the effect on the outflows of varying the 
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Figure 1. Angle-integrated mass outflow rate in unbound mate- 
rial at r = 10 9 cm (solid lines) and net mass accretion rate at the 
ISCO (dashed lines) as a function of time for model S_def (black). 
Also shown are models that suppress the nuclear binding energy 
of alpha particles (P-rec, red), self-irradiation (P-irr, green), or 
both (P-rec-irr, blue). The orbital time at r = -Ro is 2.9 ms. 



torus mass (S-m0.01, S-m0.10), torus radius (S-r75), black 
hole mass (S-M10), initial electron fraction (S-y0.05 and S- 
y0.150), initial entropy (S-s6 and S-slO), magnitude of the 
viscous stress tensor (S-v0.01 and S-0.10), and resolution 
(S-hr). Models in this series are evolved for 1000 orbits at 
r — Ro, except S-def, which is evolved for 3000 orbits at the 
same location. 

The P-series explores the effect of suppressing phys- 
ical processes relative to model S-def, to study the effect 
on the outflow generation. In model P-rec, we remove the 
binding energy of alpha particles by setting Q a — in equa- 
tion (|12[) in both the initial condition and the subsequent 
evolution. The series also explores the effect of shutting 
down self-irradiation (P-irr), as well as the effect of suppress- 
ing both recombination and self-irradiation simultaneously 
(P-rec-irr). We also test the effect of the optical depth cor- 
rections to self-irradiation, by setting Ti,i rr = (eq. [24]) in 
model P-tau (the global suppression of neutrino emission by 
Tpcak [eq. [23] is still included). All models in this series are 
evolved for 3000 orbits at r — Ro- 
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Figure 2. Evolution of the angle-averaged density in model S_def, 
illustrating the generation of the outflow. Curves correspond to 
times chosen in multiples of the orbital time at r = Ro (2.9 ms), 
showing 0, 3, 10, 30, 100, 300, 600, 1000, 1300, 1600 orbits, as 
labeled. The bold curve segments correspond to material that 
has positive energy (eq. |29|). The horizontal dotted line marks 
the density at which source terms are suppressed for numerical 
reasons, and the dashed line shows the equilibrium torus before 
relaxation for the initial condition 0)3.40 . 
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4 RESULTS 

4.1 Outflow generation 

All tori undergo the same basic evolutionary stages. The 
disk spreads in radius due to angular momentum transport, 
and the internal energy is modified by viscous heating, nu- 
clear recombination, and neutrino source terms. The accre- 
tion rate at the ISCO reaches a peak magnitude at time 
~ 0.1/q orbits at r = Ro, later decaying as a power- law in 
time after ~ 1/a orbits at the same location (Figure [T]). 

The inner regions of the disk are initially sufficiently hot 
for neutrino cooling to balance viscous heating. In contrast, 
the outer regions are colder, and viscous heating dominates 
over neutrino cooling (e.g., MPQ09). This net heating, com- 
bined with the transport of angular momentum to larger 
radii, causes an expansion in the outer regions of the disk as 



Figure 3. Ratio of neutrino heating to viscous heating Q v isc = 
T : T/(up 2 ) (eq. [TO]) in model S-def-hr at time 0.058s (orbit 20 
at r = Ro)- The white contour shows the surface Qnet = 0, inside 
which neutrino cooling dominates. 



shown in Figure [2] for model S-def. The outermost regions 
achieve positive specific energy 



- 1 2 , 1 g 

etot- 2 v P + 2(rsm )2 



+ e int + $ 



(29) 



and are therefore unbound from the central BH. 

The mass outflow rate in unbound material at r = 
10 9 cm for model S-def is shown in Figure[l] It rises from low 
values to a peak at ~ 1 s or (~ 400 orbits at r = Ro) that 
exceeds the magnitude of the mass loss at the ISCO. This 
outflow is much denser than that generated by the initial 
expansion of the equilibrium torus in the low-density ambi- 
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Figure 4. Color maps of density (left), electron fraction (middle), and entropy (right) in model S-def-hr at time 1.16s (orbit 400 at 
r = Ro), illustrating the wind morphology. Upper panels show zoomed in regions of the panels directly below, as indicated by the boxes. 
An animated version of this figure is available in the online journal. 



ent medium, and which generates the background power-law 
density profile (Figure [2|. 

The relative contribution of viscous heating, nuclear re- 
combination, and self-irradiation to the outflow generation 
is examined in Figure [T] The mass accretion rate and outflow 
rate of the default model S-def is compared with those from 
models that suppress the nuclear recombination energy (P- 
rec), self- irradiation (P-irr), and both nuclear energy and 
self-irradiation simultaneously (P-rec-irr). The recombina- 
tion of alpha particles has the strongest effect on the disk 
evolution, by (1) causing a more energetic and massive out- 
flow to arise earlier, and (2) by strongly suppressing the 
mass accretion rate at the ISCO. This drop in the accretion 
rate a t late times is consistent with the results of iLee et alj 
l|2009h . 

Self-irradiation has only a minor effect on the energet- 
ics of the system. In model S-def, neutrino heating is always 
smaller than viscous heating by at least a factor of several. 
This sub-dominance of self-irradiation can be traced to a few 



causes: (1) the neutrino luminosities are high only for a lim- 
ited period of time, because both the mass and the temper- 
ature of the torus decreases with time; (2) self-irradiation is 
the strongest right at the time where the torus becomes op- 
tically thin, and affects mostly vertical, low-density regions 
of the inner disk; and (3) the temperature in the emission re- 
gion decreases with time, resulting in lower mean energies of 
the neutrinos and thus a decreasing energy deposition rate. 
For illustration, a snapshot of the ratio of neutrino heating 
to viscous heating in model S-def-hr is shown in Figure [3] 
The effect is not entirely negligible however, as turning off 
optical depth corrections leads to ~ 20% changes in the elec- 
tron fraction and the asymptotic velocity f M4.2[) . 

The total mass accreted through the ISCO radius and 
the mass with e to t > ejected beyond a radius of 10 9 cm 
are given in Table [l] as fractions of the initial torus mass. In 
most models, approximately 10 per cent of the torus mass is 
ejected, with the rest accreting onto the BH. One exception 
is the model S-r75 with a larger initial torus radius, which 
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ejects twice as much mass as the default model. Another ex- 
ception is model S-M10, which while having a large disk ra- 
dius, also has a higher BH mass, which leads to smaller mass 
ejection by ~ 50%. The other two outliers are the models 
that suppress nuclear recombination (P-rec and P-rec-irr), 
for which ~ 3 per cent of the mass is ejected (c.f. Figure [T]). 
Doubling the resolution in each coordinate direction leads 
to insignificant quantitative differences, so we consider our 
results converged. Given the range of disk masses explored, 
the ejecta mass from delayed outflows can thus range from 
10" 3 to fO" 2 M . 

The morphology of the ejecta is quasi-spherical, with 
most of the material being ejected within ~ 60° of the disk 
midplane, as shown in Figure [4] The bulk of the disk has 
moderate entropy (~ 20), particularly in regions relevant 
for nucleosynthesis f i]4.2[) . At large radii (~ fO 9 cm), the 
outflowing material develops shocks that raise its entropy. 
Large scale instabilities are also seeded by convection in the 
disk, which cannot cool efficiently once its accretion rate 
falls below Mi gn (eq. [5]). The funnel within a few tens of 
degrees from the polar axis has different properties; we do 
not study them here given our very approximate treatment 
of general relativity, and our omission of energy deposition 
by neutrino pair annihilation. 



4.2 Outflow properties 

At early times and in the innermost regions of the disk 
(r < fOO km), weak interactions are faster than the viscous 
time, and the electron fraction reaches its equi librium value 
given the local density and temperature (e.g.. iBeloborodovl 
2003). As the disk evolves, the temperature and density de- 
crease to the point that weak interactions operate slower 
than the viscous time and the electron fraction freezes out, 
becoming an advected quantity only (MPQ09). This process 
is illustrated for our fiducial model in Figure [5] 

In the inner regions of the disk, the equilibrium elec- 
tron fraction is initially small because electrons are more 
abundant than positrons when the disk is degenerate. How- 
ever, with time F e rises to a maximum value ~ 0.3 because 
the disk becomes less degenerate as the density decreases. 
Most of this high-Y e material is accreted onto the BH, how- 
ever. The outer parts of the disk that eventually become 
unbound never achieve such high electron fractions and are 
instead ejected with Y e < 0.2. 

At a given radius, the wind material is not homogeneous 
in composition, but instead shows some variation around 
mean values. Figure [6] shows mass-flux weighted quantities 
crossing a spherical surface at the fixed radius -Rout = 1 9 cm 
in the models S-def and P-irr. Quantities are computed ac- 
cording to 



(A(R out ,t)) 



■in 



sin 9A6 Rl ut pv r A(R ou 



,*) 



f e max sinfldl? R% ut pv r 



(30) 



where A(R ou t,8 ,t) is a generic scalar quantity, and 
[#min,#max] is the angular range of the meridional integral, 
chosen as a wedge within 60° of the midplane (c.f. Figure [4}. 
The bulk of the ejected material has Y e ~ 0.2 and velocities 
v r ~ 20,000 kms" 1 . 

H eavy nuclei generatio n occurs in regions where T ~ 5 x 
10 9 K JHoffman et al.lll997T l. In the fiducial model S-def, this 




log ]0 r [cm] 



Figure 5. Angle-averaged, mass-weighted electron fraction at se- 
lected times for model S-def, as labeled. The dark-filled region 
corresponds to the weak equilibration time i WGa k = Y e /|r ne t| be- 
ing smaller than the curve time, while the light-filled region is 
such that t wca k is shorter than the viscous time t v isc- 




time [s] 

Figure 6. Instantaneous mass-flux-weighted quantities at a ra- 
dius Rout = 10 9 cm for models S-def, P-irr (no self-irradiation), 
and P-tau (no optical depth suppression). Shown are the density 
p (a), electron fraction Y e (b), and velocity v r (c). The average is 
conducted over polar angles within 60° of the equator. 
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Figure 7. Thermodynamic properties of disk material for model S-def at a radius ii nuc = 400 km, where T ~ 5 X 10 K. The histograms 
are constructed by considering all material that crosses this radius, within 60° of the midplane, over the entire disk evolution. Shown 
are the electron fraction Y e (a), entropy S (b), expansion time t exp = r/v r (c), and final helium mass fraction (d) computed from 
equation JUT}. 



corresponds to a range of radii 200 — 400 km on the midplane 
depending on the phase in the disk evolution. To quantify 
the thermodynamic properties relevant for nucleosynthesis 
in the disk, we compute time- and angle-integrated mass- 
weighted mean values 



A(R nuc ) 



JdtJ sin Ode Rl uc pv r A(R nuc , e, t) 
J dt J sin 0d6 R^ uc pv r 



(31) 



where the radius R nU c is chosen so that T ~ 5 x 10 K. The 
resulting mean values of the electron fraction, entropy, and 
expansion time at R nU c (within 60° of the midplane) are 
shown for all models in Table [l] 

The composition is largely insensitive to variations in 
the model parameters, with typical mean values Y e ~ 0.2, 
S ~ 20A:b per baryon, and expansion timescales t oxp = 
Rnuc/v r ~ 0.1 s. The magnitude of the spread from the 
mean value is illustrated in Figure [7] which breaks up the 
mean values of Y e , entropy, and expansion time for model 
S-def. 

The value of Y e is relatively insensitive to the initial 
electron fraction because most of the disk mass achieves 
weak equilibrium, erasing its initial composition. The largest 
difference in Y e occurs due to changes in the initial entropy 
of the disk and the magnitude of the viscosity. The entropy 
of the torus is related to the level of degeneracy, which in 
turn dete rmines the value of the equilibrium electron frac- 
tion (e.g.. lBeloborodovN2003l ). The magnitude of the viscos- 
ity determines the ratio of the weak equilibration time i wea k 
to the viscous times t v isc- A smaller viscosity implies that 
these two timescales become equal at a large radius, which 
increases the amount of processed Y e that enters the outflow. 
This modest dependence of Y e on a is, however, mitigated 
by two facts: (1) the disk is regulated by neutrino cooling 
to moderate degeneracy (CB07), independent of the mag- 
nitude of viscous heating; (2) the disk becomes radiatively 
inefficient (and hence thicker and less degenerate) on the 
same timescale that Y e freezes out. This coincidence results 
from the f act that e~/e + captu res both cool the disk and 
change Y e JMetzger et alj|2009eft . 

The effects of self-irradiation on the ejecta composition 



are illustrated by Figure [6] which compares the polar mass- 
flux weighted composition for model S-def, P-irr (which 
suppresses self- irradiation), and P-tau (which turns off op- 
tical depth corrections to self-irradiation). Including self- 
irradiation with optical depth correction leads to a nearly 
identical outcome than suppressing irradiation altogether. 
As with the small energetic relevance of neutrino heating, 
changes in the composition due to neutrino irradiation are 
important only around the time when the disk becomes 
transparent to neutrinos and its luminosity is still high. This 
timescale of peak irradiation is somewhat shorter than that 
required to significantly change the composition by weak in- 
teractions. The largest possible effect of neutrino heating on 
the outflow is probed by model P-tau, which completely sup- 
presses optical depth corrections to the self-irradiation lu- 
minosity (while still globally suppressing the cooling rates). 
The outflow begins earlier, achieves asymptotic velocities 
that are ~ 20% larger, and has a mean Y e that is 0.05 larger 
than the default model. Such moderate changes given an 
unrealistically large neutrino absorption rate reaffirms our 
conclusion that self-irradiation is only a small correction to 
the overall dynamics for low mass disks. 

The expansion velocity at larger radii is remarkably con- 
stant, sensitive only to the inclusion of the nuclear binding 
energy of alpha particles. 



5 DISCUSSION 

5.1 Nucleosynthesis in Disk Outflows 

Appendix[C] reviews the conditions for heavy element nucle- 
osynthesis in hot outflows, such as those from NS merger ac- 
cretion disks, and provides analytic formulae for the fraction 
of the ejected mass synthesized into heavy r-process nuclei 
Xh and the fraction X a — 1 — Ah that remains in a-particles 
(eq. [CI] ). Given the low electron fraction Y e < 0.2 of the 
disk outflows found in our calculations, we conclude that 
most of the ejecta will go into 2nd and 3rd peak r-process 
elements with mass number A ~ 130—200 (Ah ~ 1). A much 
smaller fraction, Xh c ~ 10~ 2 , will remain in 4 He, which may 
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have implications for the spectroscopic signatures of these 
events fC2)l, 

The average Galactic production rate of r-process nuclei 
due to accretion disk outflows from binary NS mergers is 
given by 



M r 



TZ NS 2 f ei M t 
lO-'Meyr" 1 



^■NS 2 \ ( /ej 

10- 4 yr-! ) \0.1 



M t 



io- 2 m q 



(32) 



where 7?. NS 2 is the rate of NS-NS mergers scaled to its esti- 
mated value i n the Milky Way (unce rtain by at least an order 
of magnitude; Kalogera et al.ll2004r ). M t is the average mass 
of the initial accretion torus, and / e j ~ 0.1 is the fraction of 
M t ejected in outflows, scaled to the characteristic value de- 
rived from our numerical simulations. Equation ((32} should 
be compared with the 'observed' rate M° bs ~ 5 x 10 _7 M o 
yr _1 required to explain the abundances of heavy A > 130 r- 
proce ss elements produced over the age of the Galaxy (jQianl 
[2000J). This shows that NS merger disk outflows consistitute 
a potentially significant r-process source. 

Although the idea that NS mergers are a promisin g 
r-process source is not new (jLattimer fe Schramm! ll974T ) , 
previous studies have focused almost exclusively on nucle- 
osynthesis of the dynamical ejecta. The disk outflows dis- 
cussed here are, however, in some ways more robust. Al- 
though the formation of a remnant torus is an agreed upon 
outcome of most modern simulations of NS-NS mergers, dis- 
agreement remains regarding the quantity— a nd even the ex- 
isten ce of— dynamically unbound material l|Faber fc Rasid 
I2012T) . 

Our simulations show no clear evidence for outflows 
powered primarily by neutrino heating. Such a wind was pre- 
viously anticipat ed to dominate mass loss from small radii 
in the disk (e.g. JMetzger et al.l I2OO8I ; ISurman et all l20Qgt 
IWanaio fc Jankal 120120 . Because neutrinos and antineutri- 
nos from the disk have similar luminosities and mean ener- 
gies, neutrino absorption is expected to drive Y e to a value 
~ 0.5. Although our simulations do show unbound polar 
outflows, the electron fraction of this material is increased 
only slightly by neutrino irradiation (Fig. [6]), since viscous 
heating dominates the unbinding of matter from the disk 
(Fig. [3]). This hierarchy is contingent upon our assumption 
of an a— viscosity, which may not reflect the true vertical 
distribution of turbulent dissipation in the disk. Our approx- 
imate treatment of neutrino physics prevents us from defini- 
tively ruling out a strong role of self-irradiation. However, if 
an even larger fraction of the dissipation occu rs in the disk 
coron a in more realistic magnetized disks (e.g. iHirose et al.l 
|2006T ). then the dominance of 'viscous'-driven (low Y e ) polar 
outflows may turn out to be a robust feature of the launching 
mechanism. Regardless of the composition of the [possibly 
neutrino-driven] winds from small radii in the disk, the late 
outflows powered by a recombination and viscous heating 
almost certainly dominate the total (time-integrated) mass 
loss. 



5.2 Radioactively-Powered Emission 

Radioactive decay of r-process elements synthesized in the 
disk outflows gives rise to an electromagnetic transient simi- 



lar to a dim supernova fe.g. lLi fc Paczvh ski 1998). Although 
the disk outflows are mildly anisotropic (Fig. [4}, the ejecta 
will become increasingly spherical as it expands homolo- 
gously with a characteristic velocity v r ~ 2 x 10 4 km s _1 
(Table [l}. Most of the energy released by the r-process oc- 
curs on a timescale ~ seconds, but this heating is lost to adi- 
abatic expansion since the outflow is highly optically thick 
at this early stage. Nevertheless, radioactive decay continues 
to add energy to the ejecta, at an approximately constant 
rate per logarithmic time (e.g. iMetzger et al.ll2010l ). Electro- 
magnetic emission peaks only once the timescale for photons 
to diffuse through the ejecta id oc nM c j/cr becomes shorter 
than the e xpansion timescale £ e x D = r/v r . This occurs on a 
timescale JMetzger et alj|20ld eqs. [3], [4]) 



6d 



t'r 



2 x lO-^kms- 1 
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1/2 



X 



10cm 2 g~ 
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(33) 



resulting in emission with a characteristic peak luminosity 



3 x 10 40 ergs _1 



M C : 



2 x lO-ikms" 1 
1/2 



-1/2 



10- 2 M P 



10cm 2 g _1 



-1/2 



(34) 



where k is the opacity, scaled to a value characteristic of 
the line opacities of Lanthanide elemen ts, given their ex- 
pected abundance in t he r-process ejecta IIKasen et alJl2013l : 
iBarnes fc Kaserj|2013h . 

Because the composition of the disk outflows (predom- 
inantly heavy r-process nuclei) is similar to that of the dy- 
namically ejected tidal tails, it may be observationally chal- 
lenging to distinguish their contribution to the transient 
emission. One unique feature of the disk outflows, however, 
is the presence of a modest amount of helium (XHe ~ 10 -2 ), 
which is not present in the dynamical ejecta due to its much 
lower entropy S < k b baryon -1 (since Xu_ c oc S 3 ' 2 ; eq. [CI]). 

Even a small quantity of He may be detectable if emis- 
sion/absorption lines (e.g. He I 1.083/im and He I 2.058/^m) 
are produced as the result non-thermal excitation, induced 
by energetic 7-rays, /3-decay electrons and fission products 
from decaying r-process elements. This is in analogy to 
the mechanism responsible for producing He lines in core- 
collapse SNe as the result of mi xing radioactive 56 Ni into 
the o uter layers of the ejecta (e.g. ILucvII1991| . iDessart et al.l 
I2OI2J) . Though a potentially promising diagnostic of condi- 
tions close to the merger site, future work is required to 
quantify the strength of He emission lines in kilonovae. 



5.3 Effects of a Long-Lived Hypermassive NS 

In the case of NS-NS mergers, our calculations implicitly 
assume that a BH forms within a few dynamical times af- 
ter coalescence. If, however, the collapse is significantly de- 
layed due to, e. g., thermal pressure and /or differential rota- 
tional support (Morrison et al.ll2004l: [O'Connor fc Ottll201ll : 
iPaschalidis et al.1 12012 ; iLehner et al.ll2012f) . then the effects 
of neutrino irradiation from the hyper-massive NS (HMNS) 
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could qualitatively alter the properties of the disk evolution 
and its outflows. 

To significantly alter the disk dynamics, neutrino emis- 
sion from the HMNS would need to maintain both a high 
flux level, and a high mean neutrino energy, for a timescale 
longer than both (1) the time required to achieve neu- 
trino transparency; and (2) the thermal time of the disk, 
such that the internal energy is raised sufficiently to un- 
bind the disk. Such neutrino-driven outflows necessarily oc- 
cur via thermal evaporation (e.g. iMc Laughlin & Surman 
20051; iMetzger et alll2008l ; ISurman et al.ll2008l ; lDessart et alj 



2009) , as in the case of proto-neutron star winds, since mo- 
mentum transfer from neutrinos is negligible [fj 

Even if neutrino irradiation from the HMNS is insuffi- 
cient to alter the dynamics of the disk, it may still alter its 
composition, especially in outer parts of the disk where the 
rate of electron/positron captures is slow. To the extent that 
the HMNSs formed in binary NS merger are similar to the 
proto-neutron stars formed in core collapse supernovae, the 
neutrino and antineutrino luminosities and temperatures are 
expected to be similar. The net effect of HMNS irradiation 
on the disk would thus be to raise Y" e in both the disk m id- 
plane and in subsequ ent outflows (IDessart et al.l 120091 ; see 
IMetzger et~aH l2009bl and iDarbha et al.ll2010l for an analo- 
gous scenario following the accretion-induced collapse of a 
white dwarf) . Spectroscopic or photometric signatures of el- 
ements synthesized exclusively in proton-rich ejecta of kilo- 
novae (e.g., Fe or 56 Ni) could thus be taken as evidence for 
a long-lived HMNS. One significant difference between an 
r-process-rich and an Fe-rich outflow is the opacity of the 
ejecta, which drastically cha nges the rise time and peak lu- 
minosity of the light curve (JBarnes fe Kasenll2013l ). Future 
work is necessary to quantify how long such a HMNS must 
survive in order to appreciably alter the outflow properties 
from the purely BH accretion scenario studied in this paper. 



5.4 Comparison with Previous Work 

iRuffert fe Jankal (119991 ) first studied the evolution of disks 
formed after NS coalescence in 3D using a realistic equa- 
tion of state, a finite volume hydrodynamic method, New- 
tonian and pseu do-Newtonian poten tials, and a neutrino 
leakage scheme ilRuffert et al.1 ll996T). Foll o w up work by 
ISetiawan et alj (J200J) and ISetiawan et al.l (|200£f ) explored 
the effect of including an a viscosity over a wide range of 
disk masses, finding that neutrino emission is sensitive to the 
magnitude of a. Our results cannot be directly compared to 
theirs, however, because their models were evolved only up 
to 70 ms. Also, due to the equation of state used, their den- 
sity floor was set at 10 8 g cm -3 , which would make it hard 
to develop any viscously driven outflow (c.f. Figure |2J. 

The interplay between magnetoh ydrodynamic and gen - 
eral relativistic effects was explored bv lShibata et al.l (J2007T) , 
finding that if poloidal fields and differential rotation are 
present, the magnetorotational instability drives turbulence 
that transports angular momentum and heats the disk. Due 
to the MHD algorithm employed, these models also require 



a floor of density ~ 10 g cm , which again prevents the 



5 Even a neutrino luminosity as large as L„ 
sub-Eddington. 



10 53 erg s x is 



development of viscously driven winds out to large radius. 
T he effect of gray radiation tr ansport has been investigated 
bv lShibata fe Sekiguchil (|2012l ). and general relativistic sim- 
ulations of t ori formed in BH-NS mergers have recently been 
rep orted by Deaton et alj (J2013J). 

iLee et al.1 (2005) evolved tori over timescales up to a few 
100 ms, using initial conditions from a merger simulation. 
They used a smoothed particle hydrodynamics method, and 
included a realistic equation of state, neutrino optical depth 
corrections, and the effect of nuclear recombination in the 
energetics of the disk. In contrast to our models, they as- 
sumed that the electron fraction satisfies beta equilibrium, 
pa rameterized the NSE ab undances by the analytic relation 
of iQian fe Wooslevl i|1996f ). included all components of the 
viscous st r ess te nsor, and focused on relatively massive tori. 
ILee et al.1 (J2009J) extended the evolution to a few seconds, 
considered a pseudo-Newtonian potential, and explored a 
wide range of tori masses and viscosity parameters. They 
found that the accretion rate drops sharply at late times 
as a consequence of including the recombination energy of 
alpha particles, an effect that we reproduce in our models 
(c.f. Figure [TJ. 

The properties of the disk evolution and outflows found 
here are in qualitative agreement with the ID model of 
MPQ09. MPQ09 predicted a somewhat larger mass of un- 
bound material, and a distribution of electron fraction ex- 
tending to higher values Y e > 0.3. In hindsight, this discrep- 
ancy can be understood by the fact that MPQ09 counted all 
mass in the disk outside the point of Y e freeze-out as eventu- 
ally becoming unbound, whereas our calculations show that 
only the outer edge of the disk (where Y e is lowest) achieves 
sufficiently high energy to become unbound. The temper- 
ature and density of the disk midplane in our simulations 
are also somewhat lower than in the ID a-disk models of 
MPQ09 (they find an inner disk temperature of ~ 10 MeV), 
resulting in less effective neutrino cooling. This discrepancy 
can partially be understood by the fact that a-disk models 
assume that viscous heating exactly balances cooling, which 
is not everywhere achieved in our models, leading the ID 
models to overestimate the midplane density. Since radia- 
tion pressure is moderately important, lower density implies 
a lower temperature in vertical hydrostatic balance. 



6 CONCLUSIONS 

We have explored the properties of outflows generated on 
the viscous time by accretion disks formed in NS-NS and 
NS-BH mergers, with an eye on detectable electromagnetic 
counterparts of the gravitational wave emission and re- 
process nucleosynthesis. Two-dimensional, time-dependent 
hydrodynamic simulations with a realistic equation of 
state, weak interactions, self-irradiation in moderately 
optically-thick environments, and angular momentum 
transport via an a-viscosity are used to characterize the 
driving, composition, and parameter dependencies of these 
outflows. Our main findings are as follows: 

1. - Disks eject ~ 10% of its mass in an unbound outflow. 
The composition is neutron rich, with typical electron 
fractions Y e ~ 0.2. This value arises because the ejected 
material comes from the outer part of the disk, which is 
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never processed to a high equilibrium electron fraction 
by weak interactions. The high- Ye material in the inner 
disk regions is accreted onto the BH (Fig. 0). Neutron-rich 
freeze-out is robust because (a) neutrino cooling regulates 
the disk to be mildly degenerate, which in turn regulates 
the equilibrium Y e prior to the freeze-out; (b) outflows 
begin once the disk becomes advective and a— particles 
form, both which occur nearly simultaneous with freeze-out. 



fe.g. lMetzger et al.ll2010f ) on the disk evolution and outflows, 
should also be explored. More accurate neutrino physics and 
three-dimensional modeling would be needed if precise pre- 
dictions for the outflow properties are necessary. 

A more detailed analysis of the nucleosynthesis prod- 
ucts and electromagnetic signal generated by these outflows 
will be carried out in subsequent work. 



2. - The outflow composition is relatively robust 
relative to initial disk parameters (Table [T]). In regions 
where heavy elements form (T ~ 5 x 10 9 K), we obtain 
characteristic entropies ~ 20 kb baryon -1 and expansion 
times ~ 0.1 s. The composition is not homogeneous, with 
spatial as well as temporal variations around the mean 
value (Figs. [4] [7]). The expected nucleosynthetic output 
consists in 2nd/3rd peak r-process elements (A > 130), 
with a small trace of He (~ 1% by mass, but > 10% by 
number), which may generate a spectroscopic signature. 

3. - Small variations ~ 10% in the electron fraction are 
introduced by changes in the initial torus entropy and the 
magnitude of the viscosity. The entropy of the outflow is 
moderately sensitive to the initial torus mass and its initial 
entropy. The amount of mass ejected is sensitive only to the 
radius at which the torus mass resides initially. The outflow 
velocity is remarkably robust, with a characteristic value of 
~ 20, 000 km s" 1 at 10 9 cm. 

4. - In addition to viscous heating, nuclear recombina- 
tion is a fundamental component in the generation of the 
outflow (Fig. [TJ. Its exclusion results in much lower ejected 
masses and asymptotic velocities. We also witnes s the drop 
in the accretion rate at late times first seen by iLee et al.l 

mi). 

5. - Self-irradiation appears to be both energetically 
and compositionally sub-dominant. This stems in part from 
the fact that irradiation is maximal around the time when 
the disk becomes transparent to neutrinos, which for most 
of our low-mass disks corresponds to early times in the disk 
evolution. The subsequent drop in disk temperatures and 
densities lowers both the amount of neutrinos emitted and 
their mean energies, diminishing the relative importance 
of neutrino energy deposition as time passes. The effects 
of self-irradiation become noticeable only when completely 
removing optical depth corrections (Fig. [6}. Even then, 
changes in global quantities are only of the order of 10%. 

There are many possible improvements to the model. 
Self-consistent angular momentum transport by magnetohy- 
drodynamic turbulence, magnetocentrifugal winds, and/or 
self-gravity (for more massive disks) would yield a first- 
principles estimate on the spatial distribution of angular mo- 
mentum and heating in the disk. A full general-relativistic 
code would help address whether the properties of the in- 
ner disk are important, particularly regarding neutrino self- 
irradiation with rapidly spinning black holes. More realistic 
initial conditions would shed light on whether the initial 
angular momentum distribution is of consequence for the 
timing and asymptotic properties of the ejecta. The possi- 
ble effects of late 'fall-back' accretion of tidal tail material 
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APPENDIX A: ACCRETION RATE FOR 
OPAQUENESS TO NEUTRINOS 

Here we provide an analytic derivation of the critical ac- 
cretion rate M opaque above which the inner disk becomes 
opaque to neutrinos. A steady-state disk is optically thick 
at radii interior to the point at which the vertical optical 
depth to neutrino absorption obeys 



t„ = £re„/2 > 1) 



(Al) 



where £ is the disk surface density and k„ is the mean opac- 
ity. For an accretion disk in steady-state, the accretion rate 
and surface density are related by M = SttvYi at radii much 
greater than the inner boundary, where v = ac 3 H is the 
kinematic viscosity; c s is the sound speed; H — c s /Qk is 
the disk scaleheight; and Q,k = (GMbh/V 3 ) 1//2 is the Kep- 
lerian orbital frequency. At accretion rates M > Mi gn , the 
inner disk is neutrino-cooled, such that the midplane is suf- 
ficiently dense that gas pressure exceeds radiation pressure, 
and the opacity is dominated by the absorption of neutrinos 
on free nuclei (CB07). In this case c s ~ (4fc s T/3m p ) 1/2 
and k„ ~ K T 2 (where kq ~ 4.5 x 10~ 39 K~ 2 cm 2 g _1 ; 
e.g. |Pi Matteo et al1l2002r i. Equation (|A1|I can thus be ex- 
pressed as a condition on the accretion rate 



M> 



87rafc_B 



() 7 ( A2 ) 

Now, for an optically thin disk, the midplane tempera- 
ture is determined by the balance between viscous heating 
gvisc = (9/4) vQ, K and neutrino cooling q„ ~ q„,oT (due to 
e~/e + captures on free nuclei, where g„,o — 8 x 10 _43 T 6 erg 
s-V'K" 6 ), viz. 



l.25(ak B /m p ) n K ' q vQ 
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. 9x io«K(^y /5 ( 
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, lY 
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-1/5 



(A3) 

-3/10 

(A4) 



where we have now expressed radius in terms of gravitational 
radii _R g = GMbh/c 2 . Substituting equation ()A4[) into equa- 
tion (IA2I). we find 
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Figure Bl. Snapshots of the angle-integrated neutrino emissivity 
for model S-def. Times shown correspond to orbits 1 (black), 3 
(red), and 10 (blue) at r = Rq. The vertical dashed line shows 
the location of i? em (eq. |16p in each of the models. 




emitting 
ring 
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Figure B2. Geometry for self-irradiation. All of the neutrino 
luminosity is assumed to originate in a ring on the midplane with 
radius i? cm . A point in the disk with coordinates (r, 6) lies at 
a distance D from an emitting element E with azimuthal angle 
</>'. The total emission received by the point in the disk is a sum 
over the emission from each ring element, leading to an angular 



dependence iang (eq. IB15I , Figure lB3l 
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Now, since most of the total neutrino luminosity is released 
from radii just outside the innermost stable circular orbit 
^isco, we can evalulate equation (|A5[) at r ~ 2Ri BCO to define 
a critical accretion rate 
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taneous volume integral of the emissivities, and with an en- 
ergy spectrum following a Fermi-Dirac distribution with zero 
chemical potential. The number density of neutrino species 
i per unit energy e, propagation direction k, and angular 
position on the ring (/>' is 

d 3 n vt 



9/5 

"""■ ! (A6) 



above which the majority of the radiated energy will orig- 
inate from an optically-thick environment. This expression 
agrees reasonably well with CB07, who find that for a — 0.1, 
Mopaquo = 0.7M Q s" 1 and 0.06M Q s" 1 for a BH of spin 



= 0(i? isco = 6R g ) and a = 0.95(R is 



1.5-R S ), 



spectively, although the dependence on a that they find, 
M p aq u oc a is slightly different than than the oc a 4 ' 5 de- 
pendence that we find in equation (]A6[) . 



APPENDIX B: SELF- IRRADIATION AND 
WEAK INTERACTIONS 

The implementation of charged-current weak inte r action 
and self-irradiation follows closely that of lFernandea (|2012l ) 
in the context of core-collapse supernovae, with suitable 
modifications to account for the disk geometry. We consider 
only electron-type neutrinos and a ntineutr i nos, w ith emis- 
sion and absorption rates following iBruennl (|1985l ). 

Our prescription for self-irradiation assumes that neu- 
trinos are emitted from a ring of material at a radius R B m- 
This approximation is motivated by the fact that the neu- 
trino emissivity displays a single peak inside the disk, where 
most of the emission is generated. Figure IBT1 illustrates the 
shape of the function for model S-def. We choose the emis- 
sion radius to be an average weighted by the neutrino emis- 
sivity (eq. [16]). This location coincides within ~ 10% with 
the radius inside which 50% of the neutrino emission is gen- 
erated. 

Each fluid element on the ring is assumed to emit 
isotropically, with a total luminosity equal to the instan- 
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where 9k is the angle between the propagation vector and 
the radial direction, and O is the step function. The intensty 
is assumed to be azimuthally symmetric around the propa- 
gation direction. The Fermi-Dirac function is given by 



Fp D (e,r, fi) = 



,(e- M )/fcT + ]_' 



(B2) 



and the minimum angle between k and the radial direction 

is 
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to lowest order in 1Z/D, where 1Z is the radius of a spherical 
emitting element in the ring, and 
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is the distance to a point in the disk with coordinates (r, 9) 
from the emitting blob in the ring, as depicted in Figure [B2l 
Keeping the lowest order terms arising from equation (|B3|I 
ensures that TL scales out of the problem. By using equa- 
tion (|B4p to capture the global geometry of the problem, we 
are ignoring light bending and Doppler effects. 

The rates per baryon of production and destruction of 
electrons and positrons that enter equation (|14f) are given 
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where j Vi and K 1 ^ i are the emissivity and absorption coeffi- 
cient, respectively, associated with electron- type neutrinos 
or antineutrinos (as subscripted). The coefficients are given 
by (|Bruennlll985h 
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where (x e is the chemical potential of electrons, n n and 
n p the number density of free neutrons and protons, re- 
spectively, Gf = Gf/(!k) 3 the Fermi constant, <?v and <?a 
the vector and axial coupling constants, respectively, and 
A m = (m n — m p )c 2 the difference between the rest mass 
energy of neutrons and protons. 

The energy source terms that enter equation (|15|) are 
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The anisotropy of the global radiation field decouples 
from the energy integrals, and enters through a dimension- 
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Figure B3. Angular dependence of the radiation field (eq. IB150 , 
The angular factor approaches a pure inverse square law at large 
distances from the origin (eq. IB16I ). 
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in all terms containing an f V( factor (D is given by eq. |B4| ). 
The spatial dependence of 7 an g is shown in Figure IB3I In 
the inner regions of the disk, this function varies by fac- 
tors of order unity following the annular geometry of the 
emission region. It becomes increasingly spherical for larger 
radii. Quantitatively, this can be seen from equation (|B16|) . 
which is a series expansion in R ern /r. Due to the azimuthal 
averaging, the leading error term scales like the 4th power 
of this ratio. The integral diverges when r — R am . To avoid 
excessive heating by artificial concentration of the emission 
at this radius, and given that the spatial distribution of the 
emission has a finite width (Figure [Blj) , we impose J ang ^ 1. 
The effect of this upper limit can also be seen in Figure [B3l 

All terms containing a f Vi factor are also multiplied by 
a normalization constant N v - 



N Vi 



L» t 



(7/16)47ri? c 2 m(TSB r4' 



(B17) 



which ensures that the incident flux is consistent with the 
total neutrino luminosity L Vi emitted by the disk. The de- 
nominator normalizes the assumed Fermi-Dirac spectrum in 
f Vi . The neutrino temperature used is the same for both neu- 
trino species, and is equal to the average, mass- weighted gas 
temperature in an (angle integrated) annular region centered 
on .Rom, with a radial width of 10%. This averaging is per- 
formed to smooth over stochastic fluctuations in localized 
disk regions. 

Weak rates are tabulated as a function of the gas tem- 
perature T, electron degeneracy parameter including rest 
mass fie/iksT), and neutrino temperature T v . The angu- 
lar function 7 ang is also tabulated in cylindrical coordinates, 
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covering the range (r/7? om ) 6 [0, 10]. Outside of this range, 
the asymptotic expansion in equation (IB16|I is used. 



APPENDIX C: HEAVY ELEMENT 
NUCLEOSYNTHESIS 

Here we review the conditions for nucleosynthesis of heavy 
elements in hot outflows, such as those from binary NS 
merger accretion disks. Free nuclei recombine into q- 
particles once the temperature decreases to T < 10 10 K, 
the energy from which plays a key role in powering the disk 
outflows explored in this paper. Heavier elements start to 
form once the temperature decreases further, T < 5 x 10 9 
K, such that the reaction 4 He(an,7) 9 Be(a,n) 12 C occurs (at 
higher temperatures, the abundance of 9 Be is limited by 
photodisintegration) . 

Immediately after 12 C forms, multiple additional a- 
captures produce heavy 'seed' nuclei with characteristic 
mass A ~ 90 — 120 and charge Z ~ 35 (the 'a-process'; 
IWooslev fc Hoffman! [l992j). Whether nucleosynthesis pro- 
ceeds to even heavier r-process nuclei depends primarily on 
the ratio of free neutrons to seed nuclei once the a-process 
completes. Since the formation of 12 C is the rate limiting 
step in forming seeds, the quantity and distribution of r- 
nuclei nuclei synthesized depends on the entropy S, electron 
fraction Y e , and expansion timescale i exp at the radius where 
T ~ 5 x 10 9 K. 

If the electron fraction of the outflow is less than that 
of the seed nuclei themelves, i.e. Y e < -j ~ 0.3, then the 
He(an,7) J Be(a,n) C reaction is ultimately limited by the 
total number of a-particles (the neutron abundance can be 
considered fixed for purposes of calculating the 12 C reaction 
yield). In this case the fina l mass fraction of alp ha particles 
is approximately given by l|Hoffman et al.lll997l ) 



Y < ~^- 



A Hc = [{2Y e y 2 +T] 
where 
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where T is proportional to the reaction rate for He — > 12 C 
integrated over the timescale available for burning ~ t C x P . 
The reaction rate depends on S 1-3 since p oc 1/S at fixed 
temperature, and since 4 He(an,7) 9 Be(a,n) 12 C is an effec- 
tive four-body reaction. Any mass not trapped in a-particles 
ends up in r-process nuclei with mass fraction 



X T = 1 — Xue, 

and mass number greater than 



A = A\\ 



2Y e 



Ji. r 



(C3) 



(C4) 



If F <C 1 then very little C forms ('He-rich freeze- 
out'), in which case helium achieves its maximum abundance 
Ahc = 2Ye, unchanged from its value just after a recombi- 
nation. Since most of the available protons are trapped in 
4 He, the r-process nuclei produced in this case are very mas- 
sive, easily reaching beyond the third peak A > 195 (and 
likely undergoing fission). Even if T = oo (X r — 1), how- 
ever, nucleosynthesis will still reach the Lanthanide elements 
(A > 139 = y4.La) if the outflow is sufficiently neutron-rich: 



0.30 — 0.40. Thus, given the low electron frac- 
tion of the disk outflows that we find in i]4.2l (Y e ~ 0.2), we 
expect L athanides to easily form, resulting in a high optical 



expect ijatnamdes to easily : 
opacity IjKasen et alj|2013i ). 
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